Data inladen

KMIdata<-read.table(file = "~/GitHub/vespa_analyses/Input/KMIdata.txt", header=TRUE, sep="\t")

# Outlier vliegsnelheid v 10m/s weglaten
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

KMIdata<-KMIdata%>%  filter(ID!=195)

library(ggplot2)
library(ggpubr)

Volgende datasets zijn gecreëerd:

KMIdata: dataset met alle vliegtijden en parameters (KMI)

Voor het model: Flight time ~ Distance gebruiken we de dataset vliegtijden per individu Omdat we hiermee de theoretische regel 1min=100m kunnen verifiëren. De imkers nemen hiervoor altijd kortste meting

Voor modellen met weerparameters, gewicht, urbanisatie: Hiervoor nemen we telkens de hele dataset, omdat elke meting van deze factoren afhangt.

We gaan verder met de dataset zonder de outliers in Melle

1.1) ForagingSpeed vs Temperature

Graphs

Van temperatuur ook categorische variabele maken (aanraden van Prof.Vangestel)

KMIdata$Temperature_cat<-cut(KMIdata$Temperature, c(8,13,17,21,25,29,33))

KMIdata$Temperature_cat

##   [1] (21,25] (21,25] (21,25] (21,25] (21,25] (21,25] (21,25] (21,25] (21,25]
##  [10] (21,25] (21,25] (21,25] (21,25] (21,25] (21,25] (21,25] (21,25] (17,21]
##  [19] (21,25] (17,21] (17,21] (17,21] (21,25] (21,25] (25,29] (25,29] <NA>   
##  [28] <NA>    <NA>    <NA>    <NA>    <NA>    (17,21] (25,29] (21,25] (17,21]
##  [37] (17,21] (25,29] (25,29] (25,29] (25,29] (25,29] (25,29] (25,29] (29,33]
##  [46] (25,29] (25,29] (25,29] (21,25] (21,25] (8,13]  (8,13]  (8,13]  (8,13] 
##  [55] (17,21] (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (17,21] (17,21]
##  [64] (17,21] (17,21] <NA>    (17,21] (17,21] (17,21] (17,21] (17,21] (17,21]
##  [73] (17,21] (17,21] (17,21] (17,21] (17,21] (17,21] (17,21] (17,21] (17,21]
##  [82] (25,29] (25,29] (25,29] (21,25] (21,25] (13,17] (13,17] (21,25] (21,25]
##  [91] (21,25] (21,25] (21,25] <NA>    <NA>    (17,21] (21,25] (21,25] <NA>   
## [100] <NA>    <NA>    <NA>    <NA>    <NA>    (13,17] (13,17] (13,17] (8,13] 
## [109] (8,13]  (8,13]  (25,29] (25,29] (13,17] (13,17] (17,21] (13,17] (13,17]
## [118] (13,17] (13,17] (13,17] (17,21] (17,21] (13,17] (13,17] (8,13]  (8,13] 
## [127] (8,13]  (13,17] (13,17] (17,21] (17,21] (13,17] (17,21] (17,21] <NA>   
## [136] (17,21] (13,17] (13,17] (13,17] (8,13]  (8,13]  (13,17] (13,17] (13,17]
## [145] (17,21] (17,21] (13,17] (17,21] (13,17] (13,17] (8,13]  (8,13]  (13,17]
## [154] (8,13]  (8,13]  (8,13]  (8,13]  (21,25] (21,25] (21,25] (17,21] (13,17]
## [163] (13,17] (13,17] (13,17] (13,17] (13,17] (13,17] (13,17] (13,17] (13,17]
## [172] (17,21] (17,21] (17,21] (17,21] (13,17] (17,21] (13,17] (13,17] (17,21]
## [181] (13,17] (8,13]  (13,17] (13,17] (13,17] (17,21] (17,21] (13,17] (13,17]
## [190] (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13] 
## [199] (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13] 
## [208] (17,21] (17,21] (17,21] (13,17] (13,17] (17,21] (17,21] (13,17] (13,17]
## [217] (17,21] (8,13]  (8,13]  (8,13]  (13,17] (13,17] (13,17] (13,17] (13,17]
## [226] (13,17] (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13] 
## [235] (8,13]  (8,13]  (13,17] (13,17] (13,17] (13,17] (13,17] (8,13]  (8,13] 
## [244] (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13]  (8,13] 
## [253] (13,17] (13,17] (17,21] (8,13]  (13,17] (13,17] (17,21] (17,21]
## Levels: (8,13] (13,17] (17,21] (21,25] (25,29] (29,33]

ggplot(data=subset(KMIdata, !is.na(Temperature_cat)), aes(x=Temperature_cat)) + geom_bar(fill="lightcoral")

Model Output

Temperature = continuous variable

Significant
Normality niet ok!

Temperature = discrete variable

Een aantal klassen significant, interpretatie?
Normality niet ok!

Poisson verdeling ook geprobeerd met glmer maar foutmelding (zie email)

library(lmerTest)

## Loading required package: lme4

## Loading required package: Matrix

## 
## Attaching package: 'lmerTest'

## The following object is masked from 'package:lme4':
## 
##     lmer

## The following object is masked from 'package:stats':
## 
##     step

library(lme4)

model_temp_cont<-lmer(Flighttime_min ~ Temperature_KMI + (1|Individualcode), offset=Distance, na.action=na.omit, data=KMIdata)
summary(model_temp_cont)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ Temperature_KMI + (1 | Individualcode)
##    Data: KMIdata
##  Offset: Distance
## 
## REML criterion at convergence: 1898.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8374 -0.2430 -0.0059  0.1753  3.8318 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 85390.122 292.216 
##  Residual                       3.288   1.813 
## Number of obs: 230, groups:  Individualcode, 91
## 
## Fixed effects:
##                  Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)     -416.0575    30.7203   91.0285 -13.543  < 2e-16 ***
## Temperature_KMI    0.3677     0.1250  138.0950   2.942  0.00382 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Temprtr_KMI -0.075

res<-residuals(model_temp_cont)
shapiro.test(res)

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.83032, p-value = 3.923e-15

qqnorm(res)
qqline(res)

model_temp_cat<-lmer(Flighttime_min ~ Temperature_cat + (1|Individualcode), offset=Distance, data=KMIdata)
summary(model_temp_cat)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ Temperature_cat + (1 | Individualcode)
##    Data: KMIdata
##  Offset: Distance
## 
## REML criterion at convergence: 1864.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9276 -0.2865 -0.0056  0.1522  3.8469 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 86380.491 293.906 
##  Residual                       3.094   1.759 
## Number of obs: 230, groups:  Individualcode, 91
## 
## Fixed effects:
##                         Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)            -412.7019    30.9886   89.0891 -13.318  < 2e-16 ***
## Temperature_cat(13,17]    2.8778     0.6471  135.0087   4.447  1.8e-05 ***
## Temperature_cat(17,21]    2.6153     0.9093  135.0332   2.876  0.00468 ** 
## Temperature_cat(21,25]    2.8244     1.5162  135.0699   1.863  0.06466 .  
## Temperature_cat(25,29]    2.4914     1.6347  135.0732   1.524  0.12983    
## Temperature_cat(29,33]  102.5019   295.5400   89.0012   0.347  0.72954    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) T_(13, T_(17, T_(21, T_(25,
## Tmp_(13,17] -0.014                            
## Tmp_(17,21] -0.020  0.469                     
## Tmp_(21,25] -0.019  0.282  0.600              
## Tmp_(25,29] -0.018  0.261  0.556  0.734       
## Tmp_(29,33] -0.105  0.001  0.002  0.002  0.002

anova(model_temp_cat, ddf="Satterthwaite", type=3)

## Type III Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)   
## Temperature_cat 64.298   12.86     5 122.3  4.1557 0.001603 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

res<-residuals(model_temp_cat)
shapiro.test(res)

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.83917, p-value = 1.037e-14

qqnorm(res)
qqline(res)

1.2) Flight time vs Distance + Temperature

Graph

Klopt dit? Kan dit ook met mixed model?

library(knitr)
library(kableExtra)

## Warning in !is.null(rmarkdown::metadata$output) && rmarkdown::metadata$output
## %in% : 'length(x) = 2 > 1' in coercion to 'logical(1)'

## 
## Attaching package: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

library("plot3D")

x <- KMIdata$Distance
y <- KMIdata$Temperature_KMI
z <- KMIdata$Flighttime_min

fit <- lm(z ~ x + y, na.action=na.exclude)
x.pred <- seq(min(x[!is.na(x)]), max(x[!is.na(x)]), length.out = 20)
y.pred <- seq(min(y[!is.na(y)]), max(y[!is.na(y)]), length.out = 20)
xy <- expand.grid( x = x.pred, y = y.pred)
z.pred <- matrix(predict(fit, newdata = xy), 
                 nrow = 20, ncol = 20)
fitpoints <- predict(fit)

scatter3D(x, y, z, pch = 19, cex = 0.6, colvar=FALSE, col="dodgerblue3", theta = 210, phi = 10, bty="u", col.panel ="grey93", expand =0.4, col.grid = "white", xlab = "Distance", ylab = "Temperature", zlab = "Flight time", surf = list(x = x.pred, y = y.pred, z = z.pred,  
facets = TRUE, col=ramp.col(col = c("dodgerblue4", "seagreen2"), n = 100, alpha=0.8), fit = fitpoints, border="black"),main = "Flight time vs Distance + Temperature")

Met plotly

Had met deze link problemen met expand.grid https://stackoverflow.com/questions/38331198/add-regression-plane-to-3d-scatter-plot-in-plotly Dan maar grid dat van hierboven gebruikt, maar nu klopt er precies iets niet? Alle datapunten liggen boven het oppervlak.

library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

library(tidyverse)

## ── Attaching packages
## ───────────────────────────────────────
## tidyverse 1.3.2 ──

## ✔ tibble  3.1.8     ✔ purrr   1.0.1
## ✔ tidyr   1.2.1     ✔ stringr 1.5.0
## ✔ readr   2.1.4     ✔ forcats 1.0.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand()          masks Matrix::expand()
## ✖ plotly::filter()         masks dplyr::filter(), stats::filter()
## ✖ kableExtra::group_rows() masks dplyr::group_rows()
## ✖ dplyr::lag()             masks stats::lag()
## ✖ tidyr::pack()            masks Matrix::pack()
## ✖ tidyr::unpack()          masks Matrix::unpack()

KMIdataNoNA<-KMIdata %>%  filter(Flighttime_min!="NA")
KMIdataNoNA<-KMIdataNoNA %>%  filter(Temperature_KMI!="NA")

fig <- plot_ly(KMIdataNoNA, x = ~Distance, y = ~Temperature_KMI, z = ~Flighttime_min, size=1)
fig <- fig %>% add_markers()
fig <- fig %>% layout(scene = list(xaxis = list(title = 'Distance'),
                     yaxis = list(title = 'Temperature_KMI'),
                     zaxis = list(title = 'Flight time (min)')))
p2 <- add_trace(p = fig,
                z = z.pred,
                x = seq(2100, 0, by = -100),
                y = seq(0, 30, by = 10),
                type = "surface")
p2

Model Ouput

Normality niet ok!
Niet significant

model_tempdist_all<-lmer(Flighttime_min ~ Distance + Temperature_KMI + (1|NestID) + (1|NestID:BaitID) + (1|Individualcode), na.action=na.exclude, data=KMIdata)

## boundary (singular) fit: see help('isSingular')

summary(model_tempdist_all)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ Distance + Temperature_KMI + (1 | NestID) +  
##     (1 | NestID:BaitID) + (1 | Individualcode)
##    Data: KMIdata
## 
## REML criterion at convergence: 1005.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2425 -0.4262 -0.1822  0.1635  4.4765 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Individualcode (Intercept) 0.1443   0.3799  
##  NestID:BaitID  (Intercept) 1.8910   1.3751  
##  NestID         (Intercept) 0.0000   0.0000  
##  Residual                   3.2814   1.8115  
## Number of obs: 230, groups:  Individualcode, 91; NestID:BaitID, 68; NestID, 34
## 
## Fixed effects:
##                   Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)      3.5399788  0.9074782 66.9650223   3.901 0.000225 ***
## Distance         0.0056150  0.0008062 62.3323406   6.965 2.42e-09 ***
## Temperature_KMI -0.0686672  0.0445457 91.2883267  -1.541 0.126655    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Distnc
## Distance    -0.336       
## Temprtr_KMI -0.903 -0.010
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

anova(model_tempdist_all, ddf="Satterthwaite", type=3)

## Type III Analysis of Variance Table with Satterthwaite's method
##                  Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Distance        159.173 159.173     1 62.332 48.5081 2.415e-09 ***
## Temperature_KMI   7.797   7.797     1 91.288  2.3762    0.1267    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

res<-residuals(model_tempdist_all)
shapiro.test(res)

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.7916, p-value < 2.2e-16

qqnorm(res)
qqline(res)

2) ForagingSpeed vs Cloudcoverage

Graphs

ggplot(KMIdata, aes(x=Cloudcoverage_KMI, y=ForagingSpeed)) + geom_point(color="lightblue") + geom_smooth(method="lm", formula =y ~ x, se=TRUE, fullrange=FALSE, level=0.95, color="darkblue") + ggtitle("All data KMI")+ theme(plot.title = element_text(hjust = 0.5, size=10))

Model Ouput

Licht significant
Normality niet ok!

model_cloud_all<-lmer(Flighttime_min ~ Cloudcoverage_KMI + (1|Individualcode), offset=Distance, na.action=na.omit, data=KMIdata)
summary(model_cloud_all)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ Cloudcoverage_KMI + (1 | Individualcode)
##    Data: KMIdata
##  Offset: Distance
## 
## REML criterion at convergence: 1574.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8876 -0.2653 -0.0075  0.1399  3.4796 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 79726.616 282.359 
##  Residual                       4.034   2.009 
## Number of obs: 191, groups:  Individualcode, 74
## 
## Fixed effects:
##                   Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)       -393.232     32.826   73.019 -11.979   <2e-16 ***
## Cloudcoverage_KMI    1.296      0.965  116.030   1.343    0.182    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Cldcvrg_KMI -0.011

anova(model_cloud_all, ddf="Satterthwaite", type=3)

## Type III Analysis of Variance Table with Satterthwaite's method
##                   Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## Cloudcoverage_KMI 7.2713  7.2713     1 116.03  1.8024  0.182

res<-residuals(model_cloud_all)
shapiro.test(res)

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.82333, p-value = 5.935e-14

qqnorm(res)
qqline(res)

3.1) ForagingSpeed vs Windspeed

Graphs

Telkens voor de modellen:

ForagingSpeed ~ Windspeed

ForagingSpeed ~ Windspeed²

plot19<-ggplot(KMIdata, aes(x=WindSpeed_KMI, y=ForagingSpeed)) + geom_point(color="khaki3") + geom_smooth(method="lm", formula =y ~ x, se=TRUE, fullrange=FALSE, level=0.95, color="olivedrab") + ggtitle("ForagingSpeed ~ Windspeed  KMI")+ theme(plot.title = element_text(hjust = 0.5, size=10))

plot22<-ggplot(KMIdata, aes(x=WindSpeed_KMI, y=ForagingSpeed)) + geom_point(color="khaki3") + geom_smooth(method="lm", formula =y ~ I(x^2), se=TRUE, fullrange=FALSE, level=0.95, color="olivedrab") + ggtitle("ForagingSpeed ~ Windspeed²  KMI")+ theme(plot.title = element_text(hjust = 0.5, size=10))


ggarrange(plot19, plot22 + rremove("x.text"), 
          labels = c("A", "B"),
          ncol = 2, nrow = 1)

Model Output

Beide modellen significant
Normality voor beiden niet ok!

model_wind_all<-lmer(Flighttime_min ~ WindSpeed_KMI + (1|Individualcode), offset=Distance, na.action=na.omit, data=KMIdata)
summary(model_wind_all)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ WindSpeed_KMI + (1 | Individualcode)
##    Data: KMIdata
##  Offset: Distance
## 
## REML criterion at convergence: 1891.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8936 -0.2627 -0.0058  0.1435  3.4706 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 85306.430 292.073 
##  Residual                       3.144   1.773 
## Number of obs: 230, groups:  Individualcode, 91
## 
## Fixed effects:
##                Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   -411.5415    30.6235   90.0651 -13.439  < 2e-16 ***
## WindSpeed_KMI    0.6332     0.1610  138.0169   3.932 0.000133 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## WindSpd_KMI -0.019

anova(model_wind_all, ddf="Satterthwaite", type=3)

## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## WindSpeed_KMI 48.614  48.614     1 138.02  15.464 0.0001325 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model_wind2_all<-lmer(Flighttime_min ~ I(WindSpeed_KMI^2) + (1|Individualcode), offset=Distance, na.action=na.omit, data=KMIdata)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues

summary(model_wind2_all)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ I(WindSpeed_KMI^2) + (1 | Individualcode)
##    Data: KMIdata
##  Offset: Distance
## 
## REML criterion at convergence: 1902.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1142 -0.2555 -0.0059  0.1360  3.6182 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 85336.122 292.123 
##  Residual                       3.302   1.817 
## Number of obs: 230, groups:  Individualcode, 91
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)        -410.12865   30.62483   90.01791 -13.392  < 2e-16 ***
## I(WindSpeed_KMI^2)    0.05677    0.01997  138.01433   2.844  0.00514 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## I(WS_KMI^2) -0.010
## optimizer (nloptwrap) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

anova(model_wind2_all, ddf="Satterthwaite", type=3)

## Type III Analysis of Variance Table with Satterthwaite's method
##                    Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)   
## I(WindSpeed_KMI^2) 26.697  26.697     1 138.01  8.0858 0.00514 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

res<-residuals(model_wind_all)
shapiro.test(res)

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.81594, p-value = 8.702e-16

qqnorm(res)
qqline(res)

res<-residuals(model_wind2_all)
shapiro.test(res)

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.80097, p-value < 2.2e-16

qqnorm(res)
qqline(res)

3.2) ForagingSpeed vs Windspeed | afhankelijk van de windrichting

Hiervoor werd telkens voor elke meting nagegaan of de hoornaar met meewind (tailwind), tegenwind (upwind) of loodrechte wind (perpendicular) te maken had. Dit volgens de formules:

|𝜃flight −𝜃wind| ≤ 45 is tailwind

45 < |𝜃flight −𝜃wind|<135 is (quasi) perpendicular

|𝜃 flight −𝜃wind| ≥135 upwind

Graphs

TO DO: Lege vakjes Wind_flight NA maken

ggplot(data=subset(KMIdata, !is.na(Wind_flight)), aes(x= Wind_flight, col=Wind_flight, y=ForagingSpeed)) + geom_boxplot()

## Warning: Removed 26 rows containing non-finite values (`stat_boxplot()`).

ggplot(data=subset(KMIdata, !is.na(Wind_flight)), aes(x= WindSpeed_KMI, col=Wind_flight, y=ForagingSpeed)) + geom_point() + facet_grid(~Wind_flight) + geom_smooth(method="lm", formula = y ~ I(x^2), se=TRUE, fullrange=FALSE, level=0.95)

## Warning: Removed 30 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 30 rows containing missing values (`geom_point()`).

Model Output

Anova van Foragingspeed en de 3 windinvloeden

interpretatie?

windanova<-lmer(Flighttime_min ~ Wind_flight + (1|Individualcode), offset=Distance, na.action=na.exclude, data=KMIdata)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues

summary(windanova)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ Wind_flight + (1 | Individualcode)
##    Data: KMIdata
##  Offset: Distance
## 
## REML criterion at convergence: 1932.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8297 -0.2558 -0.0060  0.1122  4.1052 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 84703.912 291.039 
##  Residual                       3.441   1.855 
## Number of obs: 234, groups:  Individualcode, 94
## 
## Fixed effects:
##                          Estimate Std. Error      df t value Pr(>|t|)  
## (Intercept)               -337.11     168.03   91.99  -2.006   0.0478 *
## Wind_flightperpendicular   -72.48     170.78   91.99  -0.424   0.6723  
## Wind_flighttailwind        -71.34     170.78   92.00  -0.418   0.6771  
## Wind_flightupwind          -72.06     170.78   92.00  -0.422   0.6741  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Wnd_flghtpr Wnd_flghtt
## Wnd_flghtpr -0.984                       
## Wnd_flghttl -0.984  1.000                
## Wnd_flghtpw -0.984  1.000       1.000    
## optimizer (nloptwrap) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

ForagingSpeed vs WindSpeed voor Tailwind dataset

Significant
Normality niet ok!

data_tailwind<-subset(KMIdata, Wind_flight == "tailwind")

model_tailwind<-lmer(Flighttime_min ~ WindSpeed_KMI + (1|Individualcode), offset=Distance, na.action=na.omit, data=data_tailwind)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues

summary(model_tailwind)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ WindSpeed_KMI + (1 | Individualcode)
##    Data: data_tailwind
##  Offset: Distance
## 
## REML criterion at convergence: 372.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.13557 -0.15923 -0.00299  0.01975  2.13987 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 89183.261 298.64  
##  Residual                       1.323   1.15  
## Number of obs: 42, groups:  Individualcode, 22
## 
## Fixed effects:
##                Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   -354.2999    63.6953   21.0337  -5.562  1.6e-05 ***
## WindSpeed_KMI    1.9695     0.5308   19.0056   3.710  0.00148 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## WindSpd_KMI -0.028
## optimizer (nloptwrap) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

anova(model_tailwind, ddf="Satterthwaite", type=3)

## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## WindSpeed_KMI 18.206  18.206     1 19.006  13.765 0.001484 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

res<-residuals(model_tailwind)
shapiro.test(res)

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.80132, p-value = 4.725e-06

qqnorm(res)
qqline(res)

ForagingSpeed vs WindSpeed voor Perpendicular dataset

Significant
Normality niet ok!

data_perpendicular<-subset(KMIdata, Wind_flight == "perpendicular")


model_perpendicular<-lmer(Flighttime_min ~ WindSpeed_KMI + (1|Individualcode), offset=Distance, na.action=na.omit, data=data_perpendicular)
summary(model_perpendicular)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ WindSpeed_KMI + (1 | Individualcode)
##    Data: data_perpendicular
##  Offset: Distance
## 
## REML criterion at convergence: 1086.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6773 -0.3007 -0.0055  0.1599  3.8205 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 82549.958 287.315 
##  Residual                       2.418   1.555 
## Number of obs: 132, groups:  Individualcode, 55
## 
## Fixed effects:
##                Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   -392.4157    38.7470   54.0285 -10.128 4.33e-14 ***
## WindSpeed_KMI    0.7233     0.1691   76.0084   4.277 5.46e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## WindSpd_KMI -0.016

anova(model_perpendicular, ddf="Satterthwaite", type=3)

## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## WindSpeed_KMI  44.22   44.22     1 76.008  18.289 5.456e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

res<-residuals(model_perpendicular)
shapiro.test(res)

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.84239, p-value = 1.443e-10

qqnorm(res)
qqline(res)

ForagingSpeed vs WindSpeed voor Upwind dataset

Niet significant
Normality niet ok!

data_upwind<-subset(KMIdata, Wind_flight == "upwind")


model_upwind<-lmer(Flighttime_min ~ WindSpeed_KMI + (1|Individualcode), offset=Distance, na.action=na.omit, data=data_upwind)
summary(model_upwind)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ WindSpeed_KMI + (1 | Individualcode)
##    Data: data_upwind
##  Offset: Distance
## 
## REML criterion at convergence: 511.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9621 -0.2480 -0.0058  0.1086  2.9613 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 89543.074 299.237 
##  Residual                       4.592   2.143 
## Number of obs: 56, groups:  Individualcode, 27
## 
## Fixed effects:
##                Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   -412.9125    57.6072   26.0315  -7.168 1.29e-07 ***
## WindSpeed_KMI   -0.2379     0.4400   28.0061  -0.541    0.593    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## WindSpd_KMI -0.025

anova(model_upwind, ddf="Satterthwaite", type=3)

## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## WindSpeed_KMI 1.3425  1.3425     1 28.006  0.2924  0.593

res<-residuals(model_upwind)
shapiro.test(res)

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.80469, p-value = 3.644e-07

qqnorm(res)
qqline(res)

Model Selection

fullmodel<-lmer(Flighttime_min ~ Traject100m + WindSpeed_KMI + Temperature_KMI + Cloudcoverage_KMI + Weight_ind + Wind_flight + (1|Individualcode), offset=Distance, na.action=na.omit, data=KMIdata)
summary(fullmodel)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ Traject100m + WindSpeed_KMI + Temperature_KMI +  
##     Cloudcoverage_KMI + Weight_ind + Wind_flight + (1 | Individualcode)
##    Data: KMIdata
##  Offset: Distance
## 
## REML criterion at convergence: 410.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3917 -0.3001 -0.0556  0.2722  3.3894 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 53883.948 232.129 
##  Residual                       3.547   1.883 
## Number of obs: 71, groups:  Individualcode, 16
## 
## Fixed effects:
##                      Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)         -550.3334   308.0425   13.0598  -1.787  0.09723 . 
## Traject100m         1673.6657   550.0121   12.9981   3.043  0.00943 **
## WindSpeed_KMI         -0.1017     1.6066   50.0768  -0.063  0.94978   
## Temperature_KMI        0.5069     0.6809   50.0637   0.744  0.46009   
## Cloudcoverage_KMI     -1.2680     2.4523   50.0386  -0.517  0.60739   
## Weight_ind          -479.3634   861.3756   13.0048  -0.557  0.58731   
## Wind_flighttailwind    0.9777     1.3264   50.0094   0.737  0.46447   
## Wind_flightupwind      0.4210     1.0123   49.9981   0.416  0.67929   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trj100 WS_KMI Tm_KMI Cl_KMI Wght_n Wnd_flghtt
## Traject100m -0.352                                              
## WindSpd_KMI  0.050 -0.019                                       
## Temprtr_KMI -0.052  0.017 -0.930                                
## Cldcvrg_KMI -0.031  0.016 -0.591  0.558                         
## Weight_ind  -0.874 -0.105 -0.025  0.025  0.012                  
## Wnd_flghttl -0.003  0.006 -0.220  0.001  0.303  0.001           
## Wnd_flghtpw -0.010  0.006 -0.188  0.158  0.296  0.003  0.238

model1<-lmer(Flighttime_min ~ Traject100m + Temperature_KMI + Cloudcoverage_KMI + Weight_ind + Wind_flight + (1|Individualcode), offset=Distance, na.action=na.omit, data=KMIdata)
summary(model1)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ Traject100m + Temperature_KMI + Cloudcoverage_KMI +  
##     Weight_ind + Wind_flight + (1 | Individualcode)
##    Data: KMIdata
##  Offset: Distance
## 
## REML criterion at convergence: 413.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4009 -0.3049 -0.0583  0.2752  3.4189 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 53842.509 232.040 
##  Residual                       3.479   1.865 
## Number of obs: 71, groups:  Individualcode, 16
## 
## Fixed effects:
##                      Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)         -549.3616   307.5411   13.0053  -1.786  0.09737 . 
## Traject100m         1673.0213   549.7055   12.9996   3.043  0.00942 **
## Temperature_KMI        0.4667     0.2481   51.0125   1.881  0.06569 . 
## Cloudcoverage_KMI     -1.3589     1.9593   51.0134  -0.694  0.49111   
## Weight_ind          -480.7074   860.7830   12.9992  -0.558  0.58603   
## Wind_flighttailwind    0.9596     1.2813   51.0107   0.749  0.45734   
## Wind_flightupwind      0.4090     0.9845   51.0036   0.415  0.67961   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trj100 Tm_KMI Cl_KMI Wght_n Wnd_flghtt
## Traject100m -0.352                                       
## Temprtr_KMI -0.016 -0.001                                
## Cldcvrg_KMI -0.002  0.006  0.029                         
## Weight_ind  -0.874 -0.106  0.006 -0.004                  
## Wnd_flghttl  0.008  0.002 -0.568  0.220 -0.004           
## Wnd_flghtpw  0.000  0.003 -0.046  0.233 -0.002  0.206

model2<-lmer(Flighttime_min ~ Traject100m + Temperature_KMI + Cloudcoverage_KMI + Weight_ind + (1|Individualcode), offset=Distance, na.action=na.omit, data=KMIdata)
summary(model2)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ Traject100m + Temperature_KMI + Cloudcoverage_KMI +  
##     Weight_ind + (1 | Individualcode)
##    Data: KMIdata
##  Offset: Distance
## 
## REML criterion at convergence: 418
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5700 -0.2799 -0.0505  0.2098  3.4708 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  Individualcode (Intercept) 53760.67 231.863 
##  Residual                       3.39   1.841 
## Number of obs: 71, groups:  Individualcode, 16
## 
## Fixed effects:
##                    Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)       -551.1020   307.2951   13.0051  -1.793  0.09618 . 
## Traject100m       1671.8808   549.2846   12.9996   3.044  0.00941 **
## Temperature_KMI      0.5674     0.2009   53.0059   2.825  0.00666 **
## Cloudcoverage_KMI   -1.7809     1.8501   53.0103  -0.963  0.34013   
## Weight_ind        -477.7494   860.1200   13.0005  -0.555  0.58802   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trj100 Tm_KMI Cl_KMI
## Traject100m -0.352                     
## Temprtr_KMI -0.014  0.000              
## Cldcvrg_KMI -0.004  0.006  0.179       
## Weight_ind  -0.874 -0.106  0.004 -0.003

model3<-lmer(Flighttime_min ~ Traject100m + Temperature_KMI + Cloudcoverage_KMI + (1|Individualcode), offset=Distance, na.action=na.omit, data=KMIdata)
summary(model3)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ Traject100m + Temperature_KMI + Cloudcoverage_KMI +  
##     (1 | Individualcode)
##    Data: KMIdata
##  Offset: Distance
## 
## REML criterion at convergence: 1536.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5917 -0.2547 -0.0054  0.1572  3.2093 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 63054.182 251.106 
##  Residual                       3.757   1.938 
## Number of obs: 191, groups:  Individualcode, 74
## 
## Fixed effects:
##                    Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)       -559.5955    45.6541   72.7061 -12.257  < 2e-16 ***
## Traject100m       1018.4148   225.8033   72.0127   4.510 2.46e-05 ***
## Temperature_KMI      0.4846     0.1545  115.1085   3.138  0.00216 ** 
## Cloudcoverage_KMI    1.0641     0.9342  115.0364   1.139  0.25704    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trj100 Tm_KMI
## Traject100m -0.766              
## Temprtr_KMI -0.070  0.012       
## Cldcvrg_KMI -0.002 -0.001 -0.079

model4<-lmer(Flighttime_min ~ Traject100m + Temperature_KMI  + (1|Individualcode), offset=Distance, na.action=na.omit, data=KMIdata)
summary(model4)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Flighttime_min ~ Traject100m + Temperature_KMI + (1 | Individualcode)
##    Data: KMIdata
##  Offset: Distance
## 
## REML criterion at convergence: 1867.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8341 -0.2430 -0.0060  0.1773  3.8263 
## 
## Random effects:
##  Groups         Name        Variance  Std.Dev.
##  Individualcode (Intercept) 70771.693 266.029 
##  Residual                       3.288   1.813 
## Number of obs: 230, groups:  Individualcode, 91
## 
## Fixed effects:
##                  Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)     -559.5730    42.8324   89.6604 -13.064  < 2e-16 ***
## Traject100m     1009.9184   228.1899   89.0088   4.426 2.72e-05 ***
## Temperature_KMI    0.3728     0.1250  138.0877   2.983  0.00338 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trj100
## Traject100m -0.757       
## Temprtr_KMI -0.061  0.009

library(cAIC4)

## Loading required package: stats4

## Loading required package: nlme

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:lme4':
## 
##     lmList

## The following object is masked from 'package:dplyr':
## 
##     collapse

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:plotly':
## 
##     select

## The following object is masked from 'package:dplyr':
## 
##     select

library(AICcmodavg)

## 
## Attaching package: 'AICcmodavg'

## The following object is masked from 'package:lme4':
## 
##     checkConv

AIC(fullmodel)

## [1] 430.4662

AIC(model1)

## [1] 431.2467

AIC(model2)

## [1] 431.9586

AIC(model3)

## [1] 1548.611

AIC(model4)

## [1] 1877.757

Multicollinearity check

library(ggcorrplot)
library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:purrr':
## 
##     some

## The following object is masked from 'package:dplyr':
## 
##     recode

datacor<-KMIdata[ , c("Temperature_KMI", "Cloudcoverage_KMI", "WindSpeed_KMI", "Weight_ind", "Traject100m", "NestID")]
source("~/GitHub/vespa_analyses/Input/HighstatLibV10.R") 
corvif(datacor)

## 
## 
## Variance inflation factors
## 
##                       GVIF
## Temperature_KMI   3.538307
## Cloudcoverage_KMI 1.333698
## WindSpeed_KMI     3.811471
## Weight_ind        2.981154
## Traject100m       1.435785
## NestID            5.380406

cormat <- round(cor(datacor, use = "pairwise.complete.obs"), 2)
ggcorrplot(cormat, lab= TRUE, type = "lower", ggtheme = ggplot2::theme_gray,
   colors = c("#6D9EC1", "white", "#E46726"))

Masterthesis script,weerparameters KMI

Heleen Van Ransbeeck

2023-02-27

Data inladen

1.1) ForagingSpeed vs Temperature

Graphs

Model Output

1.2) Flight time vs Distance + Temperature

Graph

Model Ouput

2) ForagingSpeed vs Cloudcoverage

Graphs

Model Ouput

3.1) ForagingSpeed vs Windspeed

Graphs

Model Output

3.2) ForagingSpeed vs Windspeed | afhankelijk van de windrichting

Graphs

Model Output

Anova van Foragingspeed en de 3 windinvloeden

ForagingSpeed vs WindSpeed voor Tailwind dataset

ForagingSpeed vs WindSpeed voor Perpendicular dataset

ForagingSpeed vs WindSpeed voor Upwind dataset

Model Selection

Multicollinearity check